sf packagesf package to handle
spatial data?”Packages needed
library(dplyr)
library(sf)
library(tibble)
library(ggplot2)
library(remotes)
remotes::install_github('ropensci/osmdata')
library(osmdata)
library(osmextract)
Let’s start from where we left off:
bb <- getbb('Rotterdam', format_out = 'sf_polygon')
x <- opq(bbox = bb) %>%
add_osm_feature(key = 'building') %>%
osmdata_sf ()
buildings <- x$osm_polygons %>% st_transform(.,crs=28992)
buildings$build_date <- as.numeric(ifelse(buildings$start_date <1945, 1945,buildings$start_date))
ggplot(data = buildings) +
geom_sf(aes(fill = build_date, colour=build_date)) +
scale_fill_viridis_c(option = "viridis")+
scale_colour_viridis_c(option = "viridis")
sf?sf is a package which supports simple features (sf), “a standardized
way to encode spatial vector data.”. It contains a large set of
functions to achieve all the operations on vector spatial data for which
you might use traditional GIS software: change the coordinate system,
join layers, intersect or unit polygons, create buffers and centroids,
etc. cf. the sf cheatsheet.
We are going to go through some of these basics with the case study of Rotterdam buildings.
Let’s focus on really old building and imagine we’re in charge of their conservation. We want to know how much of the city would be affected by a non-construction zone of 500m around pre-1800 buildings.
Let’s select them and see where they are.
old <- 1800 # year prior to which you consider a building old
old_buildings <- filter(buildings, start_date <= old)
ggplot(data = old_buildings) + geom_sf()
As conservationists, we want to create a zone around historical buildings where building regulation will have special restrictions to preserve historical buildings.
Let’s say this zone should be 500 meters. In GIS terms, we want to
create a buffer around polygons. The corresponding function
sf is st_buffer, with 2 arguments: the
polygons around which to create buffers, and the radius of the
buffer.
distance <- 500 # in meters
buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
ggplot(data = buffer_old_buildings) + geom_sf()
Now, we have a lot of overlapping buffers. We would rather create a
unique conservation zone rather than overlapping ones in that case. So
we have to fuse the overlapping buffers into one polygon. This operation
is called union and the corresponding function is
st_union.
single_old_buffer <- st_union(buffer_old_buildings) %>%
st_cast(to = "POLYGON") %>% st_as_sf() %>%
add_column("ID"=as.factor(1:nrow(.))) %>% st_transform(.,crs=4326)
For the sake of visualisation speed, we would like to represent
buildings by a single point (for instance: their geometric centre)
rather than their actual footprint. This operation means defining their
centroid and the corresponding function is st_centroid.
sf::sf_use_s2(FALSE)
centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)
ggplot() +
geom_sf(data = single_old_buffer, aes(fill=ID)) +
geom_sf(data = centroids_old)
Now what we would like to distinguish conservation areas based on the
number of historic buildings they contain. In GIS terms, we would like
to know how many centroids each fused buffer polygon contains. This
operation means intersecting the layer of polygons with the layer of
points and the corresponding function is
st_intersection.
centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>%
add_column(n=1)
centroid_by_buffer <- centroids_buffers %>%
group_by(ID) %>%
summarise(
n = sum(n)
)
single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)
ggplot() +
geom_sf(data = single_buffer, aes(fill=n_buildings)) +
scale_fill_viridis_c(alpha = 0.8,
begin = 0.6,
end = 1,
direction = -1,
option = "B")
Let’s map this layer over the initial map of individual buildings.
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,
begin = 0.6,
end = 1,
direction = -1,
option = "B")
old <- 1939 # year prior to which you consider a building old
distance <- 100 # in meters
old_buildings <- filter(buildings, start_date <= old)
buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
single_old_buffer <- st_union(buffer_old_buildings) %>%
st_cast(to = "POLYGON") %>% st_as_sf() %>%
add_column("ID"=1:nrow(.)) %>% st_transform(.,crs=4326)
centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)
centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>%
add_column(n=1)
centroid_by_buffer <- centroids_buffers %>% group_by(ID) %>% summarise(n = sum(n))
single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,begin = 0.6,
end = 1,direction = -1, option = "B")
old <- 1939
distance <- 100
#select
old_buildings <- filter(buildings, start_date <= old)
#buffer
buffer_old_buildings <- st_buffer(old_buildings, dist = distance)
#union
single_old_buffer <- st_union(buffer_old_buildings) %>%
st_cast(to = "POLYGON") %>% st_as_sf() %>%
add_column("ID"=1:nrow(.)) %>% st_transform(.,crs=4326)
#centroids
centroids_old <- st_centroid(old_buildings) %>% st_transform(.,crs=4326)
#intersection
centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>%
add_column(n=1)
centroid_by_buffer <- centroids_buffers %>%
group_by(ID) %>%
summarise(
n = sum(n)
)
single_buffer <- single_old_buffer %>% add_column(n_buildings = centroid_by_buffer$n)
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,
begin = 0.6,
end = 1,
direction = -1,
option = "B")
Problem: there are many pre-war buildings and the buffers are large so the number of old buildings is not very meaningful. Let’s compute the density of old buildings per buffer zone.
single_buffer$area <- st_area(single_buffer, ) %>% units::set_units(., km^2)
single_buffer$old_buildings_per_km2 <- as.numeric(single_buffer$n_buildings / single_buffer$area)
ggplot() +
geom_sf(data = buildings) +
geom_sf(data = single_buffer, aes(fill=old_buildings_per_km2), colour = NA) +
scale_fill_viridis_c(alpha = 0.6,
begin = 0.6,
end = 1,
direction = -1,
option = "B")
We have seen how to create spatial buffers and centroids, how to intersect vector data and how retrieve the area of polygons.
In short:
sf package to treat geospatial datast_* functions for basic GIS operationsggplot package to map the results